r. 


AD-AQ91  247 
UNCLASSIFIED 


ARMY  ARMAMENT  RESEARCH  AND  DEVELOPMENT  COMMAND  ABERD— ETC  F/G  9/2 
USER'S  MANUAL  FOR  THE  BRL  SUBROUTINE  TO  EVALUATE  SINE.  COSINE.  — ETC(U) 
AUG  BO  J  N  WALBERT.  E  M  WlNEHOLT 

ARBRL-TR-02255  SBIE-AD-E430  514  NL 


m 


i  # 


fa. 


i 


\ 


TECHNICAL  REPORT  ARBRL-TR-Q2255 


USER'S  MANUAL  FOR  THE  BRL  SUBROUTINE 
TO  EVALUATE  SINE,  COSINE,  AND  EXPONENTIAL 
INTEGRALS  FOR  COMPLEX  ARGUMENT 


Jarafts  N.  Walbert 
Emma  M.  Wlnaholt 


August  1980 


DTI  r. 

|ELE 
NOV  5  1960 


Y 


US  ARMY  ARMAMENT  RESEARCH  AND  DEVELOPMENT 
BALLISTIC  RESEARCH  LABORATORY 

ABERDEEN  PROVING  GROUND.  MARYLAND 


Approved  for  public  roltost;  distribution  unllaltod. 


14  074 


f 


Dehttoynhis  report  whon  It  is  ao  longer  assist* 
Do  ast  mum  it  to  tho  originator. 


Secondary  distribution  of  this  report  by  originating 
or  sponsoring  activity  is  prohibited. 

Additional  eopios  of  this  roport  nay  bo  obtained 
fron  tho  National  TOchnieal  Infoxnation  Service, 

U.S.  Dopartaont  of  Coasorco,  Springfisld.  Virginia 
22151. 


Tho  findings  in  this  roport  are  not  to  bo  construed  as 
sn  official  Dopartaont  of  tho  Any  position,  unless 
so  designated  by  other  suthorisod  documents. 


UNCLASSIFIED 


SECURITY  CLARIFICATION  OF  THIS  PAGE  fWlMl  Data  g, HkwQ 


REPORT  DOCUMENTATION  PAGE 

READ  INSTRUCTIONS 

BEFORE  COMPLETING  FORM 

TECHNICAL  REPORT  ARBRL-TR-02255  ATVAO^,^4/*"] 

S.  RECIPIENT'S  CATALOG  NUMBER 

4.  TITLE  (mtd  SubtUlo) 

USER’S  MANUAL  FOR  THE  BRL  SUBROUTINE  TO  EVALUATE 

s.  tyre  or  report  a  period  covered 

Technical  Renort 

SINE,  COSINE,  AND  EXPONENTIAL  INTEGRALS  FOR 

COMPLEX  ARGUMENT 

S.  PERFORMING  ORO.  REPORT  NUMBER 

7.  AUTHORS) 

James  N.  Walbert 

Emma  M.  Wineholt 

a.  contract  or  orant  numbers 

S.  PERFORMING  ORGANIZATION  NAME  AND  ADDRESS 

U.S.  Army  Ballistic  Research  Laboratory 

ATTN:  DRDAR-BLP 

Aberdeen  Provine  Ground.  MD  21005 

to.  PROGRAM  ELEMENT.  PROJECT,  TASK 
AREA  *  WORK  UNIT  NUMBERS 

1L161102AH43 

II.  CONTROLLING  OFFICE  NAME  AND  AOORCSS 

U.S.  Army  Armament  Research  §  Development  Command 
U.S.  Army  Ballistic  Research  Laboratory 

^SlrdeeS^rovSiie  Ground.  MD  2100S 

u.  report  date 

August  1980 

IS.  NUMBER  OF  PAOES 

76 

14.  MONITORING  AGENCY  NAME  4  ADORESV"  dUtoront  from  Controlling  Ottteo) 

IS.  SECURITY  CLASS,  (of  tblo  report) 

tlNCLASSIP.TEII 

is*,  declassification/ downgrading 
SCHEDULE 

Is.  DISTRIBUTION  STATEMENT  (at  thla  Rapart) 

Approved  for  public  release,  distribution  unlimited 

17.  DISTRIBUTION  STATEMENT  (al  Ih.  afeatracl  antarad  In  Black  SO,  II  dlilmranl  from  Rapati) 

it.  supplementary  notes 

It.  KEY  WOROS  (Continue  on  from  o  Ido  It  noeoooary  m d  Identity  by  block  man  bor) 

Sine  Integral 

Cosine  Integral 

Exponential  Integral 

Gauss  Continued  Fraction 

Subroutine 

ZB.  ABSTRACT  (T-ifMmi  amaaaarm  BBFim— T  end.  Idcmiltf  ky  klcak  mmkac) 

The  formulas  used  in  the  BRL  Sine,  Cosine  and  Exponential  Integral 
subroutine  are  derived  in  sufficient  detail  to  enable  the  user  to  correctly 
implement  the  computer  code.  Estimates  of  precision,  sample  output,  and  an 
annotated  listing  of  the  computer  code  are  included  in  the  appendices. 

DO  i  j5Tn  U73  niTiMer  iMvwitMKLiTi  UNCLASSIFIED 


SECURITY  CLARIFICATION  OR  THIS  SUM  Tier 


TABLE  OF  CONTENTS 

Page 


NOTATION .  5 

I.  INTRODUCTION .  7 

II.  INPUT  AND  OUTPUT  VARIABLES .  8 

III.  METHODS  OF  COMPUTATION 

A.  Definitions  .  9 

B.  Series:  |z|  £  10 .  9 

C.  Gauss  Continued  Fraction:  10  <  |z|  <.  75 . 10 

D.  Asymptotic  Series:  |z|  >  75 . 12 

E.  Recurrence  Relation  for  the  Exponential  Integral.  ...  12 

F.  Programming  Methods  .  13 

IV.  CONCLUSIONS . 19 

V.  ACKNOWLEDGEMENTS . 19 

REFERENCES . 21 

APPENDICES 

A.  PRECISION  ESTIMATES  .  23 

B.  LISTING  OF  SUBROUTINE  SCINT  .  25 

C.  SAMPLE  OUTPUT  FROM  SUBROUTINE  SCINT  .  37 

D.  DERIVATIONS  OF  THE  COMPUTATIONAL  FORMULAS  .  57 

DISTRIBUTION  LIST . 75 


NOTATION 


This  user's  manual  is  designed  to  assist  the  mathematician  or 
programmer  using  the  BRL  Sine,  Cosine,  and  Exponential  Integral  Sub¬ 
routine.  FORTRAN  symbols  for  variables  and  arithmetic  operations  are 
used  in  the  body  of  the  report  for  consistency  with  excerpts  from  the 
coding . 

As  an  aid  to  the  reader  unfamiliar  with  standard  FORTRAN,  the 
following  symbols  are  defined: 


Symbol 

Operation 

Algebraic 

Notation 

FORTRAN 

Notation 

1. 

add 

a  +  b  ■ 

A  +  B 

2. 

- 

subtract 

a  -  b  * 

A  -  B 

3. 

* 

multiply 

a  x  b  * 

A  *  B 

4. 

/ 

divide 

a  *  b  =* 

A  /  B 

Numbers  are  written  in  specific  ways  to  define  their  type: 


1.  Integer:  2 

2.  Real:  2.  or  2.0  5 

3.  Standard  notation  2.78  xlO  :  2.78  E+05 

(double  precision)  2.78  D+05 


FKKIKLNO  PiOB  XUMC-MV  IlliCD 


I .  INTRODUCTION 


Sine  and  cosine  integrals  occur  in  applications  of  Tranter's 
method*  to  the  evaluation  of  stresses  in  thick -walled  cylinders.  The 
need  for  highly  precise  values  of  the  integrals  for  complex  argument 
led  to  the  development  of  this  BRL  subroutine.  Although  tables  of  sine, 
cosine,  and  exponential  integrals  have  been  published,  these  tables  are 
of  necessity  limited  in  scope.  Moreover,  interpolation  between  given 
values  in  such  tables  results  in  loss  of  accuracy. 

While  computer  subroutines  exist  for  the  computation  of  certain  of 
these  integrals  for  restricted  values  of  the  argument,  the  authors  know 
of  no  subroutine  which  is  valid  for  the  wide  range  of  complex  argument 
and  order  of  the  exponential  integral  or  which  has  the  degree  of 
precision  of  the  subroutine  presented  in  this  report. 

The  three  methods  used  to  compute  the  values  of  the  integrals  are 


1.  Series, 

2.  Gauss  continued  fractions,  and 

3.  Asymptotic  series. 


Each  of  these  methods  will  be  discussed  in  sufficient  detail  to  enable 
the  user  to  understand  the  subroutine.  Recourse  to  special  multiple 
precision  codes  or  integer  arithmetic  has  been  deliberately  avoided. 

The  goal  was  to  provide  a  high  degree  of  precision  through  careful 
attention  to  analytic  detail. 

The  subroutine  has  been  written  in  double -precis ion  FORTRAN  IV  and 
has  been  code-checked  on  the  CDC  7600.  Examples  run  on  the  CDC  7600 
have  agreed  to  25  significant  digits  with  tables  of  the  sine  integral 
generated  by  C-B  Ling2.  (See  Appendix  C) .  Complex  arithmetic  has  not 
been  used  in  the  computer  code;  the  annotated  listing  in  Appendix  B, 
together  with  the  analysis  presented  in  section  III,  will  serve  to 
illustrate  the  manner  in  which  real  and  imaginary  parts  of  the  integrals 
are  computed.  As  a  general  rule,  whenever  the  real  and  imaginary  parts 
of  a  number  are  stored  in  an  array  of  length  2,  the  real  part  is  in  the 
first  location  and  the  imaginary  part  is  in  the  second  location. 

- 

C.  J.  Tranter ,  Integral  Transforms  in  Mathematical  Physios,  Methuen  and 
Co.,  Ltd.,  London,  England,  1966. 

Chih-Bing  Ling,  Collected  Papers ,  Vol.  II,  Virginia  Polytechnic  Institute 
and  State  University,  Blacksburg,  VA,  1979, 
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II.  INPUT  AND  OUTPUT  VARIABLES 
The  subroutine  statement  is 

SUBROUTINE  SCINT  (X,  Y,  SI,  Cl,  EX,  NORDBR,  ICODE,  I ERR) . 

The  input  variables  are  X,  Y,  NORDER,  and  ICODE.  X  and  Y  are  double¬ 
precision  real  variables,  while  NORDER  and  ICODE  are  integer  variables. 
ICODE  describes  to  the  subroutine  the  manner  in  which  X  and  Y  are  to  be 
interpreted.  If 


ICODE  =  1,  the  complex  argument  z  is  x+iy 

* 

=  2,  the  complex  argument  z  is  x  EXP(iy) 

NORDER  is  the  order  n  of  the  exponential  integral  En(z)  to  be  computed. 

The  output  variables  are  SI,  Cl,  EX,  and  IERR.  SI,  Cl,  and  EX  are 
double-precision  real  arrays  of  length  2,  and  IERR  is  an  integer  variable 
used  as  an  error  code. 

SI(1)  *  Re  Si (z) 

SI (2)  *  Im  Si(z) 

CI(1)  *  Re  Ci(z) 

Cl  (2)  -  Im  Ci(z) 

EX(1)  *  Re  En(z),  n  -  NORDER 

EX(2)  «  Im  En(z),  n  -  NORDER 

IERR  *  0,  no  errors  occurred; 

*  1,  input  value  for  ICODE  was  not  1  or  2; 

=  2,  the  magnitude  of  z  was  less  than  l.D-48,  and  interpreted 
to  be  0. ; 

*  3,  the  argument  of  z  was  180°  degrees; 

»  4,  the  magnitude  of  z  was  negative; 

■  5,  negative  order  was  specified  for  En (z) ■ 

It  should  be  noted  that  not  all  of  the  nonzero  values  for  IERR 
indicate  fatal  errors.  If  IERR  -  1  or  4,  no  computations  are  performed. 

If  IERR  ■  2  (that  is,  if  z  -  0),  then  the  sine  integral  Si(0)  ■  0,  while 
the  cosine  integral  Ci(0)  is  undefined.  If  n  •  NORDER  ■  0  or  1,  E  (0) 
is  also  undefined.  If  n  ■  NORDER  >1,  E  (0)  ■  l/(n-l).  If  IERR  ■  3,  Ci(z) 
amd  En(z)  are  not  defined  for  z  on  the  negative  real  axis;  the  values  of 
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Si (z)  are  computed.  Finally,  if  IERR  ■  5,  E  (z)  is  not  computed,  but 
Si(z)  and  Ci(z)  are. 

Appropriate  error  messages  are  printed  to  accompany  the  nonzero 
values  of  IERR.  It  is  necessary  to  declare  the  double-precision  arrays 
SI (2),  Cl (2),  and  EX(2)  in  a  DIMENSION  statement  in  the  calling  program. 

III.  METHODS  OF  COMPUTATION 


A. 


Definitions 

The  sine  integral  is  defined 
z 


Si(z) 


I 


sin  t 
t 


dt 


by 


for  all  complex  z. 

The  cosine  integral  is 

z 

Ci(z)  «  y  +  In  z  + 


defined  by 


cos  t  -  1 
t 


dt 


(1) 


(2) 


for  |arg  z|<if,  where  y  *  .57721  ...  is  Euler's  constant  and  In  z  is  the 
complex  logarithm  of  z. 

In  z  =  log  [ z |  +  i  arg  z.  (3) 

The  exponential  integral  is  defined  by 

oo 

E"(z)  =  /  dt  W 

*i  z 

for  Re  z  >  0  and  order  n  =  0,1,2,... 

B.  Series:  |  z  |<_  10 

Series  expansions  for  Si(z)  and  CiC?)  are  given  by 
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!  '  ¥ 


Si(z) 


y  M r 

Zj  T57TT 


k  2k+l 
z 


Ef  llkz2k 

-fei  (iTp.  •  iar* 


For  the  exponential  integral,  one  has 


En(z) 


|arg  2  |  <tt 


£&(—•♦£  i\-± 

V  k-l  /  k-0 

k*n-l 


»(7) 


Since  the  infinite  geries  in  Eqs.  (5),  (6),  and  C7)  are  majorized  by 
an  infinite  series  for  e  ,  they  are  absolutely  convergent  throughout  the 
finite  complex  plane.  The  presence  of  the  term  In  z  in  Eqs.  (6)  and 
(7)  invalidates  these  equations  at  the  origin  and  along  the  negative  real 
axis.  As  |z|  increases,  the  rate  of  convergence  decreases.  Numerical 
experimentation  showed  agreement  between  the  series  and  the  continued 
fractions  (to  be  discussed  nextl  in  the  region  8  |z|  ^  12  to  in  excess 
of  28  significant  digits.  |z|  ■  IQ  was  therefore  chosen  as  the  cutoff 
point  for  utilization  of  the  series  expansion.  See  Appendix  D  for  the 
derivation  of  these  expansions. 

C.  Gauss  Continued  Fraction:  10< j z | <75 

For  the  exponential  integral,  the  continued  fraction  is  giyen  by 


En(z) 


valid  for  |arg  z|  <  ir.  (See  Appendix  D) 
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The  continued  fraction  expansion  in  Eq.  (8)  can  also  be  used  to 
evaluate  the  sine  and  cosine  integrals.  In  particular,  for  |arg  z|  < 

SiU)  =  21  [Ei(iz)  *  E1C-iz) J  +  J  .  (9) 

and 

Ci(z)  =  -  j^Ej(iz)  +  Ej(-iz)j  .  (10) 

(See  Appendix  D  for  the  derivations.)  Then  using  the  continued  fraction 
in  Eq.  (8)  to  evaluate  E^(±  iz) ,  one  obtains  the  values  of  Si(z)  and 

Ci(z)  from  Eqs.  (9)  and  (10),  respectively,  for  |arg  zj  <  j  and  10<)z|<75 

For  z  such  that  y  <|arg  z|<  it,  Si(-z)  and  Ci(-z)  are  computed,  and  then 

use  is  made  of  the  fact  that 


Si(z)  =  -  Si (~z)  (11) 


and 


Ci  (z)  =  Ci(-z)  +  iir  (12) 

For  |arg  z|  =  j  ,  a  problem  arises  in  Eqs.  (9)  and  (10),  since 

either  iz  or  -iz  may  lie  on  the  negative  real  axis,  the  branch  cut  for 

E1(z).  It  will  be  shown  in  Appendix  D  that  the  continued  fraction  expan¬ 
sion  in  Eq.  (8)  is  still  valid  when  z  <  0,  if  used  properly.  That  is, 

if  z  is  real  and  negative,  say  z  =  -  x,  x  >  Q,  then  define 


EjC-x)  =  -  Ei(x)  +  i  it  , 


C13) 


where 


Ei(x) 


(141 


P.V.  denotes  the  Cauchy  principal  value  of  the  integral,  and  the 
sign  of  the  iir  term  is  chosen  according  to  the  following  convention: 
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•_  a 


(15) 


lim  E. (-x  +  iy)  »  -  Ei(x)  -  iff 
y-K)+  1 


lim  E.  (-x  +  iy)  *  -  Ei(x)  +  in 
y-KT  L 


The  continued  fraction  expansion  in  Eq.  (8)  converges  to  -Ei(x) 
when  z  *  -  x,  x  <  0.  Suppose,  then,  that  z  ■  ix,  x  >  0.  In  this  case, 
iz  »  -  x,  and  the  continued  fraction  converges  to  -  Ei(x).  Since  z  has 
been  rotated  in  the  positive  (anti-clockwise)  sense  to  the  negative  real 
axis,  it  is  subtracted  from  the  imaginary  part  of  -  Ei(x)  (which  is  0) 
to  provide  the  value  of  Ej(iz).  Similarly,  if  z  ■  -  ix,  then  -  iz  ■  -  x 
by  rotation  in  the  negative  sense,  and  ff  is  added  to  the  imaginary  part 
of  -Ei(x)  to  obtain  the  value  of  E|(-iz).  Kith  these  conventions,  Eqs. 
(9)  and  (10)  may  be  used  to  compute  Si  and  Ci  on  the  imaginary  axis. 

The  value  |z|  »  7S  was  chosen  as  the  cutoff  point  for  use  of  the 
continued  fraction,  since  numerical  experimentation  showed  agreement 
between  the  continued  fraction  and  the  asymptotic  series  (to  be  discussed 
next)  to  in  excess  of  26  significant  digits  in  the  range  |z|  ■  73  to  |z| 

=  78  • 

D.  Asymptotic  Series:  |z|  >75 

For  the  exponential  integral,  the  asymptotic  series 


E.W-V  i  -  5.  .  "fell  -  n(n*lUn*2)  , 


is  valid  for  |arg  z]  <  j  ir  ,  provided  that  iff  is  added  or  subtracted  from 
the  value  of  the  series,  as  appropriate,  when  the  branch  cut  is  crossed 
(see  Appendix  D).  In  particular,  when  z  is  on  the  negative  real  axis, 
the  series  in  expression  (16)  is  asymptotic  to  -Ei(xl,  so  the  computation 
of  Si  and  Ci  from  this  asymptotic  series  is  accomplished  using  Eq.  Cal 
and  (10)  in  exactly  the  same  manner  as  was  described  in  section  C  for 
the  continued  fraction. 


Recurrence  Relation  for  the  Exponential  Integral 


For  n  ^  1,  it  is  shown  in  Appendix  D  that 
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for  all  z.  This  recurrence  relation  is  stable  for  increasing  h  whenever 
n  <  |z|,  and  is  stable  for  decreasing  n  whenever  n  >  |z|. 

F.  Programming  Methods 

Throughout  this  section,  reference  is  made  to  program  line  numbers 
which  can  be  found  in  Appendix  B.  All  computations  within  the  subroutine 
are  done  using  the  rectangular  coordinate  form  of  the  complex  quantities. 
Since  input  to  the  subroutine  may  be  in  rectangular  or  polar  form,  lines 
180  through  440  check  the  form  of  the  input,  and  convert  polar  input  to 
rectangular  form.  Polar  input  angles  are  modified,  if  necessary,  to 
values  greater  than  -  180°  and  less  than  or  equal  to  180°  by  adding  or 
subtracting  an  appropriate  multiple  of  360°.  It  should  be  noted  that 
input  variables  are  not  modified  by  the  subroutine;  in  lines  160  and  170, 
auxiliary  variables  are  assigned. 

Lines  500  through  630  comprise  the  series  computation  section.  The 
array  TSAVE,  of  length  2,  contains  the  real  and  imaginary  parts  of  the 
term 

z)N0RDER-l 

(NORDER-l) !  ’ 


for  use  in  the  exponential  integral  in  lines  620  and  630.  In  lines  580 
through  610,  the  terms  of  the  summation  in  the  quantity 


NORDER-1 


are  computed. 


Lines  660  through  880  make  up  the  continued  fraction  section.  The 
first  call  to  subroutine  CONTFR  (line  670)  is  to  compute  the  value  of 
En(z),  n  =  NORDER,  z  input.  Lines  680  through  760  compute  Ej(iz)  for  use 

in  Eqs.  (9)  and  (10),  as  described  in  section  III  C.  If  z  lies  in  the 
right  half-plane  (line  690),  z  is  rotated  90°  in  the  counterclockwise 
sense  by  setting  XT  =  -  Y  and  YT  =  X.  That  is,  if  z  =  X  +  iY,  then  iz  * 

-  Y  +  iX  (lines  730  and  740).  If  z  is  in  the  left  half-plane  or  on  the 
imaginary  axis,  z  is  replaced  by  -  z  before  rotation  by  ±  i  (lines  700 
through  720) .  Subsequent  to  the  computation.  Si  and  Ci  are  modified  in 
accordance  with  Eqs.  (11)  and  (12)  of  section  III  C  (lines  810  through 
880).  Having  computed  E^(iz)  in  line  750,  iz  is  converted  to  -iz  in 

lines  770  and  780,  and  E^(-iz)  is  computed  in  line  790.  If  the  input 

value  of  Z  was  imaginary  (X  *  0),  then  either  YT  (line  740)  or  YT  (line 
780)  will  be  0  with  the  corresponding  XT  being  negative.  As  described 


in  section  III  C,  the  imaginary  part  of  the  computed  Ej  is  modified  by 
±  ir,  as  appropriate,  in  line  760  or  in  line  800.  Lines  810  through  840 
apply  Eqs.  (9)  and  (10). 

The  asymptotic  series  computation  section  is  in  lines  900  through 
1130.  The  first  call  to  subroutine  ASYMP  (line  920)  computes  E„(z),  n  - 
NORDER,  z  input.  Lines  950  through  1130  implement  Eq.  (9),  (10),  (11), 
and  (12)  exactly  as  described  in  the  continued  fraction  section  above, 
except  of  course  that  subroutine  ASYMP  is  called  instead  of  subroutine 
CONTFR. 

Lines  1160  through  1510  make  up  the  error  handling  section  of  the 
subroutine.  For  details,  see  the  description  of  IERR,  section  II. 

Subroutine  SERIES  is  comprised  of  lines  1640  through  2380.  The 
arrays  PSS(100,2),  PSC(100,2),  and  PSE(200,2)  are  used  to  store  the 
computed  terms  of  the  series  for  SI,  Cl,  and  ER,  respectively.  The  first 
index  is  in  each  case  the  term  number;  the  second  index  corresponds  to 
the  real  (1)  and  imaginary  (2)  part  of  the  term.  The  terms  are  computed 
in  descending  order  and  summed  in  ascending  order,  to  minimize  round-off 
error.  For  the  range  of  application  of  the  series  representations,  the 
terms  are  monotone  decreasing  except  for  perhaps  the  first  two  or  three 
terms  when  |z|  >6,  and  hence  it  is  unnecessary  to  sum  the  series  from 
two  directions. 

For  SI  and  Cl,  one  requires  terms  of  the  form 


and 


(-nV*1 

(2n«li  (2n«l) !  ’ 

respectively.  Letting  z*X  ♦  iY,  these  terms  are  computed  as  follows: 


— T  =■  FACTR  +  i  FACTI 
m : 

(-Dn  »  sw 

m  -  EM 
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In  line  1740,  EM  is  initialized  as  1.D0,  and  in  lines  1760-1770,  PACTR 
and  FACTI  are  initialized  as  X/EM  and  Y/EM,  respectively.  SW  is  initially 
given  the  value  +  1.D0  (line  1800].  Lines  1820  through  2220  make  up  the 
computational  loop.  The  terms  for  SI  are  computed  first  (lines  1830- 
1840]  as 


SW*FACTR/EM 


and 


SW*FACTI/EM 


EM  is  incremented  by  1.D0  in  line  1940,  The  sign  of  SW  is  changed  in 
line  1950,  and  FACTR+iFACTI  is  multiplied  by  X/EM  +  i  Y/EM  (lines  1970- 
2010).  The  terms  for  Cl  are  then  computed  in  lines  2020-2030  as 


SW*FACTR/EM 


and 


SW*FACTI/EM 


EM  is  again  incremented  by  1.D0  (line  2150),  and  FACTR+iFACTI  is  multi¬ 
plied  by  X/EM+iY/EM  (lines  2170-2210).  Note  that  the  sign  of  SW  is  not 
changed  this  time,  since  the  next  term  in  the  series  representation  for 
SI  has  the  same  sign  as  the  term  just  computed  for  Cl.  The  loop  index 
increments,  and  the  next  term  in  the  series  representation  for  SI  is 
computed. 

Computation  of  the  series  representation  for  E  requires  terms  of 
the  form  n 


f  i  \  m  m 
(-1)  Z 

(m-n+1 ) (m! ) 


Again  letting  z=X  +  iY,  these  terms  are  computed  as 


m 

■^-r  =  FACTR+iFACTI 


(-1)®  =  ESW 

m  ■  EM 
n-1  =  EXFACT 

~rr  =  EXCOEF 
m-n+i 


Since  the  series  representation  for  E  contains  terms  both  even  and  odd 
powers  of  z,  two  terms  of  the  series  ?or  E„  are  confuted  each  time 
through  the  loop.  (Hence  the  dimension  PSE  (200, z)). 

As  indicated  elsewhere,  the  array  TSAVE,  dimensioned  2,  holds  the 
real  and  imaginary  parts  of  the  factor 

(n-ijl 


used  in  the  series  for  E  (z).  If  n-1,  TSAVE(1)«1 . 00  and  TSAVE(2)»0.D0, 
the  values  given  these  variables  initially  in  lines  1780-1790.  Once 
into  the  computational  loop,  so  long  as  mj*n-l,  that  is,  EM  f  EXFACT 
(lines  1850  and  2040),  the  terms  PSE  are  computed  as 


ESW*  EXCOEF*  FACTR 


and 


ESW*EXCOEF*FACTI 


The  sign  of  ESW  is  changed  after  each  such  computation,  and  FACTR  and 
FACTI  are  modified  as  described  above. 

If  m-n-1,  that  is,  EM  =  EXFACT,  the  term  in  z11-1  is  omitted  from 
the  summation,  so  that  PSE-0.D0  (lines  1880-1890  or  lines  2070-2080). 
In  this  case. 


TSAVE (1)=-ESW*FACTR 
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and 


TSAVE(2)=-ESW*FACTI 


(lines  1860-1870  or  lines  2050-2060).  The  minus  sign  occurs  because  the 
loop  computes  terms  of  the  series 

00 

(-l)kzk  (-l)k+1zk 

(k-n+l)k!  ‘  LmJ  (k-n+l)k!  ’ 

k-0 
k>*n-l 

and  TSAVE  is  required  to  have  the  value 


E 

k*0 

Mn-1 


,  ...k  k 
(-1)  z 

k! 


where  k=n-l. 


The  variable  TERM  is 


m 


and  is  computed  in  line  2130.  When  TERM  <  EPS,  the  computation  loop  is 
exited  (line  2140).  Note  that  if  N  <  NORDER,  the  loop  is  not  exited, 
since  in  this  case  TSAVE  has  not  yet~been  computed.  In  any  event,  no 
more  than  NORDER+1  terms  (NMAX)  are  computed.  If,  after  NMAX  terms, 
TERM  is  not  less  than  EPS,  a  message  is  printed  (line  2230).  In  lines 
2250  through  2340,  the  terms  are  summed.  If  N0RDER>1,  the  first  term 
in  the  series  representation  is  EXPACT,  which  is  added  to  EX(1)  in  line 
2350. 


Subroutine  C0NTFR  is  in  lines  2390  through  2680.  This  subroutine 
computes  the  2*NMAXt‘1  convergent  of  the  continued  fraction  in  Eq.  (8), 
where  NMAX-100.  The  variable  ADD,  initialized  to  NMAX  in  line  2420, 
will  take  on  the  values  1  to  NMAX,  in  descending  order.  The  variable 
ADDZ,  initialized  to  NORDER+NMAX  in  line  2430,  will  take  on  the  values 
NORDER  to  NORDER+NMAX,  in  descending  order.  W(ll  and  W(21  are  the  real 
and  imaginary  parts,  respectively,  of  the  value  of  the  convergent  at 
each  stage  of  computation.  W(l)+iW(2]  is  initialized  to  the  value 
(NORDER+NMAX) +i0. DO  in  lines  2440-2450. 
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-  . 


Lines  2460  through  2S60  comprise  the  mein  computational  loop. 
z*X+iY  is  added  to  W(l)+iW(2)  in  lines  2470  and  2480.  N(l)+iN(2)  is 

inverted  by  writing 


1 . DO/ (W (1 ) +iW (2) )  -  W(l)/R  -  iN(2)/R, 
where  R«W(1)*W(1)*W(2)*W(2) 


R  is  computed  in  line  2490,  and  in  lines  2500  and  2510,  the  new  values 
of  W(l)  and  W(2)  are  computed  by  inverting,  multiplying  by  the  current 
value  of  ADD,  and  adding  1.D0. 

Next,  ADDZ  is  decreased  by  1.D0,  and  new  values  of  NCI)  and  W(2) 
are  computed  by  inverting  and  multiplying  by  the  current  value  of  ADDZ. 
Finally,  ADD  is  decreased  by  1.D0  and  the  loop  begins  again.  After  the 
final  pass  through  the  loop,  W(l)  and  W(2)  contain  the  real  and 
imaginary  parts,  respectively,  of 

...  NORDER 

H  -  -  i 

1  +  NORDER+1 

z  +  - - 


1 

z+norder+nMax 


The  desired  convergent  is 


1 

z+W  * 


so,  in  lines  2570  through  2610,  z  is  added  to  N  and  the  sum  is  inverted. 
In  lines  2620  through  2660,  this  result  is  multiplied  by 


e 


-z 


e 


-X-iY 


■  e'Xcos  Y  -  ie’Xsin  Y. 


Subroutine  ASYMP  is  in  lines  2690  through  3120.  As  in  subroutine 
SERIES,  the  values  of  the  terms  are  saved  in  the  PSE  array,  and  are 
summed  in  descending  order  of  the  index.  As  in  subroutine  CONTFR,  only 
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the  function  E  is  computed.  The  asymptotic  expansion  requires  terms  of 
the  form  n 

(-l)k~1(n)(n+in...)(n+k-2) 


for  k=2, 3, . . . ,  NMAX. 


These  terms  are  computed  in  the  DO-LOOP  in  lines  2840  through  2950. 
The  Nth  term  is  obtained  from  the  previous  one  by  multiplying  by  - 
(NORDER+N-1)  and  dividing  by  z  (lines  2850  through  2870).  TERM  is  the 
magnitude  squared  of  each  term;  the  computation  loop  is  exited  prior  to 
the  computation  of  all  NMAX  terms  if  TERM  is  larger  than  |z|  ,  or  if 

TERM  is  less  than  l.D-44. 

The  terms  are  summed  in  reverse  order  in  lines  3000  through  3050. 

To  these  sums  is  added  z"*  (lines  3060  and  3070),  after  which  they  are 
multiplied  by  e-z  (lines  3080  and  3090). 

IV.  CONCLUSIONS 

The  accuracy  of  the  derivation  and  coding  of  the  computational 
formulas  used  in  this  subroutine  has  been  carefully  checked.  Because 
sufficiently  precise  tables  are  not  available,  the  only  means  of  checking 
the  subroutine  for  complex  z  is  through  alternate  computational  methods. 
As  has  been  noted,  agreement  between  series,  continued  fraction,  and 
asymptotic  expansion  in  their  regions  of  overlap  provides  some  degree  of 
verification.  In  certain  parts  of  the  complex  plane,  the  functions  have 
known  values  with  which  the  subroutine  may  be  compared.  For  example,  on 

the  imaginary  axis.  Re  Si(z)  =  0  and  Im  Ci(z)  =  ±  -j  ,  precisely  the 

values  given  by  the  subroutine  in  verification  runs.  Some  precision 
estimates  are  given  in  Appendix  A. 
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APPENDIX  A 


PRECISION  ESTIMATES 

A  full  discussion  of  error  analysis  and  alternate  Beans  of  computat- 
ion  of  Si,  Ci,  and  En  will  be  the  subject  of  a  future  report  by  the 
principal  investigator.  Preliminary  results  show  that  evaluation  by 
Gaussian  quadratures  agrees  with  the  values  obtained  by  the  continued 
fraction  to  27  significant  digits  in  double  precision  CDC  FORTRAN. 
Applications  of  the  recurrence  relation  (equation  17)  for  Ejjfz)  in  a 
Miller-type  algorithm  have  revealed  similar  agreement  with  the  series 
computations  and  the  asymptotic  expansions. 

Internally,  the  subroutine  determines  the  limits  of  computation  as 
follows:  the  series  computations  are  stopped  when  the  magnitude  of  the 
last  computed  term  is  less  than  l.D-144,  or  when  100  terms  have  been 
computed,  whichever  comes  first.  The  continued  fraction  computes  the 
100th  convergent.  This  value  was  found  to  provide  the  most  stable 
results  in  the  range  of  application.  The  asymptotic  expansion  is  term¬ 
inated  in  one  of  three  ways:  when  200  terms  have  been  computed,  when 
the  magnitude  of  the  last  computed  term  was  less  than  l.D-144,  or  when 
the  magnitudes  of  the  computed  terms  reach  a  minimum. 

Internal  computations  are  performed  using  the  rectangular  coordinate 
form  of  the  complex  quantities.  For  this  reason,  the  output  of  the  sub¬ 
routine  will  be  inherently  more  precise  when  the  input  is  in  rectangular 
form.  For  example,  an  input  of  X-0,  Y»A,  I CODE- 1,  will  produce  better 
results  than  an  input  of  X-A,  Y-90*,  ICODE-2,  because  of  both  the 

truncation  error  involved  in  converting  90#  to  j  radians,  and  the 

subsequent  buildup  of  errors  due  to  the  fact  that  the  computed  value  of 
X  may  not  be  identically  0. 
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APPENDIX  B 

LISTING  OF  SUBROUTINE  SCINT 
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SAMPLE  OUTPUT  FROM  SUBROUTINE  SCINT 
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APPENDIX  D 

DERIVATIONS  OF  THE  COMPUTATIONAL  FORMULAS 
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appendix  d 

DERIVATIONS  OF  THE  COMPUTATIONAL  FORMULAS 

The  series  expansions  in  Eqs.  (5)  and  (6)  are  derived  from  term- 
by-term  integration  of  the  series  expansions  of  the  integrands  in  the 
integrals  in  Eqs.  (1)  and  (2),  respectively. 


From  Eq.  (1) , 


Si  (z) 


-  f  L  V  c-nkt2k>1 
J  t  L  (2k+l) !  dt 
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®  z 

=  y  c-dk  c  t2k.t 
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Similarly,  from  Eq.  (2), 


Ci (z)  -  y  -  In  z  = 
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To  derive  the  series  expansion  in  Eq.  (7),  one  proceeds  as  follows: 
from  Eq.  (4) , 


co 


GO 


dt  . 


It 


is  known 


that 


Euler's  constant  y  is  given  by 


Thus 


provided  z  is  not  on  the  negative  real  axis  or  at  the  origin.  Then 


E.J.  Whittaker  and  G.N.  Watson,  A  Course  of  Modem  Analysis, 
Cambridge  University  Press,  London ,  19.27,  p.  246. 
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Ej(2)  -  -  Y  ♦ 


or, 

e,m  -  -  t  -  m  z 

k*l 


Starting  with  the  integral  in  Eq.  (4),  an  integration  by  parts 
yields  the  recurrence  relation  of  Eq.  (17): 


provided  n  >  1.  That  is, 


From  the  series  expansion  for  E^(z)  and  the  recurrence  relation  in 
Eq.  (17),  one  can  derive  the  series  expansion  for  En(z).  To  this  end, 
apply  the  recurrence  relation  (n-1)  times: 
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n-1 
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c-l)  (Z)  e  (-1)  Z  e  (-1)  Z  e 

(n-1)  (n-1)  (n-2)  (n-l) (n-2) (n-3) 


.  (-i)”-2,"-2.-1  ,  (-d"-1.-1  p ... 

(n-l)!  (n-1)!  El(z) 


Expanding  each  of  the  e  terms  in  a  power  series,  and  using  the  series 
expansion  for  E^z), 


En(z) 
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(-D  2 
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■n _ I _  L 

(n-l)!  ^ 
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:-ukzk 

TTi- 


Redefining  the  indices  of  summation,  and  collecting  all  terms  with  z  , 


En^z-*  (n-l)k!  T  ^ 

k*0  k-1 

Wn-1  kj<n-l 

(Equation  continued  on  next  page) 
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k-n-2 
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oo 

z 
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(n-1) ! (k-n+1) (k-n+1) ! 


From  this  last  equation,  one  sees  that  the  only  contributions  to 
powers  of  2  less  than  (n-1)  come  from  the  first  (n-1)  series.  In  particular, 
for  any  integer  £  such  that  0  <_  £  <_  (n-2) ,  there  are  (£+1)  terms  with  z4 
given  by 


(-1)  z 


1 


(n-1) £! 


(n-1) (n-2V(£-l) ! 


.  4* 


(n-l) (n-2) . 7 . (n- £- l) (0) T 


Repeated  factoring  of  the  expression  above  in  braces  yields 


(-i)V 
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(S*r) 

(Equation  continued  on  next  page) 
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En(z) 


n-l  n-1 


(n-1) ! 


(• 


In  z  + 


£0 


E 

k=0 

Mn-1 


k  k 


n+l)k! 


which  is  Eq.  (7) . 

The  continued  fraction  of  Eq.  (8)  is  derived  from  the  Gauss  continued 
fraction  by  a  method  found  in  Wall^.  The  Gauss  continued  fraction  is  given 
by 


F (a,b+l ,c+l ;z) 
F(a,b,c;z) 


1  - 


a(c-b) 

c(c+l) 

(b+l) (c-a+1) 

1  (c*l)  (c+2)  z 

(a+1) (c-b+1) 


1  - 


(c*2)  (c+3) 


where  F(a,b,c;z)  is  the  hypergeometric  f unction , 


F(a,b,c;z) 


(aMbH  k 

(c)kk!  Z 


and 


(a)k  - 


T (a+k) 
r(a) 


the  quotient  of  two  gamma  functions. 

If,  in  the  series  for  the  hypergeometric  function,  one  replaces  z 
by  -  cz  and  takes  the'  limit  as  c  -*■  the  divergent  series 


4 

H.  S.  Wall ,  Continued  Fractions.  D.  Van  Nostrand  Co.,  Inc.,  New  York, 
1948 ,  pages  336-352. 
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xt*.  #41J  . 


■v-rtmuz 


fl(a,b;-z)  »  1  -  abz  +  a(a+l) (b) Cb*l)  |y 


3 

-  a(a+l)(a+2)(b)(b+l)(b+2)  —■  +... 


is  obtained.  Using  the  same  transformation  and  limit  in  the  continued 
fraction  of  Gauss, 
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!2(  a,b;-z) 
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Using  the  divergent  series  for  fl. 
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(1+zu)  e  u  du. 


where  use  has  been  made  of  the  binomial  coefficient 


(0 


_! 


(j-k)!(k!) 


•  rvn  = 


ft') 


and  of  the  integral  representation 

oo 

F(p)  =  e_uu^"^du. 


It  follows  that 
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Choosing  b=0  and  using  the  continued  fraction  expansion  for  the  quotient 
on  the  left. 


Since  fl(a,b;-z)  =  Sl(b,a;-z) 
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Setting  b=l, 
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It  can  be  shown  that  the  integrals  in  this  last  equation  converge 
for  all  values  of  z  not  on  the  negative  real  axis.  Replace  z  by  ^  in  the 
first  of  these  integrals  and  let  t  =  1  +  —  . 
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-u 


,.  1  .a 

(!  +  r  u) 


du 


z 


ze 


or,  on  letting  a  =  n, 


This  continued  fraction  may  be  simplified  to  the  form 


for 

|arg  z\  <  ■n ,  which  is  Eq.  (8). 

To  derive  Eqs 

.  (9)  and  (10),  note 

that  for 
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Moreover, 
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The  first  term  in  square  brackets  above  tends  to  0 
5 

as  s  -*■  <»,  since 


/{* 


-  cos  t 


The  second  term  in  square  brackets  is 


♦  In  s  -  In  z  +  In  z 


w.  Magnue,  F.  Oberhettinger  and  R.P.  Soni ,  Formulae  and,  Theorems 
for  the  Special  Funotione  of  Mathematical  Physioe*  Springer  Verlag , 
New  York,  1988,  page  35. 
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4 


ts>|=l 


proving  Eqs.  (9)  and  (10). 


It  remains  to  derive  the  asymptotic  expansion  for  En(z).  To  this 
end,  rearrange  the  recurrence  relation  in  Eq.  (17)  in  the  form 


En(x) 


e'z  -  nWz)) 


and  use  it  repeatedly: 
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Therefore, 
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2  V  '  z  z2  z3  ) 


This  derivation  is  equivalent  to  repeated  integration  by  parts  starting 
with  the  integral  in  Eq.  (4).  (see,  e.g.,  Olver6,  page  67). 

Q  ' 

F.  V.  J.  Olver,  Asymptotics  and  Special  Functions,  Academic  Press,  New 
York,  1974. 
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All  of  the  computational  formulas  used  in  this  report  can  be  found 
in  reference  7.  It  should  be  noted  that  this  reference  contains  several 
other  computational  formulas  which  would,  at  first  glance,  seem  to 
provide  a  more  simple  means  of  evaluating  Si  and  Ci  for  large  |z|  than 
the  method  used  in  this  report.  In  particular,  the  sine  and  cosine 
integrals  may  be  written  in  terms  of  the  auxiliary  functions 


f(z) 


00 

/-zt 

V- dt  * 

t +1 


and 


g(z)  = 


or 

/ 


„  -zt 
~ —  dt. 
t  +1 


The  functions  f  and  g  have  asymptotic  expansions  which  are  easily  derived, 
but  which  fail  to  represent  Si  and  Ci  correctly  on  the  imaginary  axis. 

This  problem  will  be  examined  in  more  detail  in  a  subsequent  report  on 
verification  of  the  present  subroutine. 


7 

M.  Abramowitz  and  I.  Stegun,  editors.  Handbook  of  Mathematical 
Functions ,  National  Bureau  of  Standards ,  U.S.  Dept,  of  Commerce , 
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